import numpy as np  # 线性代数的一个重要的包
import pandas as pd  # 数据处理
import matplotlib.pyplot as plt  # 画图

path = "ex1data1.txt"
data = pd.read_csv(path, header=None, names=['Population', 'Profit'])
print(data.head())  # 浏览数据
# 默认输出了前五行
#    Population   Profit
# 0      6.1101  17.5920
# 1      5.5277   9.1302
# 2      8.5186  13.6620
# 3      7.0032  11.8540
# 4      5.8598   6.8233
# 其中，Population表示一个区域中有多少个人，Profit表示每个人的平均收入
# 人口对于收入的影响关系，使用线性回归来描述

print(data.describe())
# 输出更加详细的数据的描述
#        Population     Profit
# count   97.000000  97.000000      # 样本个数，共97个
# mean     8.159800   5.839135      # 平均值
# std      3.869884   5.510262      # 方差
# min      5.026900  -2.680700
# 25%      5.707700   1.986900      # 从小到大排，取第25%大的数据
# 50%      6.589400   4.562300      # 中位数
# 75%      8.578100   7.046700
# max     22.203000  24.147000


# 数据可视化
data.plot(kind='scatter', x='Population', y='Profit', figsize=(12, 8))  # kind='scatter或ling' 绘制散点图 或 直线图，figsize控制图像大小
plt.show()


# 接下来就要使用 梯度下降法 实现 线性回归

# 在数据集前插入一列，使我们可以使用向量化的解决方案来计算代价和梯度，在训练集的左侧插入一列全为'1'的列
# 以便计算x0 = 1， loc为0， name为ones，value为1
data.insert(0, 'Ones', 1)   # 在第1列，插入一列数，列名为'Ones'，值都为1
print(data.head())
#    Ones  Population   Profit
# 0     1      6.1101  17.5920
# 1     1      5.5277   9.1302
# 2     1      8.5186  13.6620
# 3     1      7.0032  11.8540
# 4     1      5.8598   6.8233


# shape[0] 表示行数，shape[1] 表示列数
cols = data.shape[1]  # 一共3列，因此 cols 的值为3
X = data.iloc[:, 0:cols - 1]  # iloc 作用：抽出某些行和某些列，在[]中用‘，’分隔 要抽取的行 和 要抽取的列
# 此处即取所有行的第一列 到 第二列
print(X)
# X如下：
#     Ones  Population
# 0      1      6.1101
# 1      1      5.5277
# 2      1      8.5186
# 3      1      7.0032
# 4      1      5.8598
# ..   ...         ...

y = data.iloc[:, cols - 1:cols]  # 取2到3列，即结果值y
print(y.head())
# y如下：
#     Profit
# 0  17.5920
# 1   9.1302
# 2  13.6620
# 3  11.8540
# 4   6.8233


# ==============================================================
# 代价函数应该是numpy矩阵，所以我们需要转换X和Y后才能使用，还需要初始化theta，即把theta所有元素都设置为0
X = np.matrix(X.values)  # 转化为矩阵
# 转化后的X如下
# [[ 1.      6.1101]
#  [ 1.      5.5277]
#  [ 1.      8.5186]
#  [ 1.      7.0032]
#   ...........
#  [ 1.      5.4369]]

y = np.matrix(y.values)
# 转化后的y如下 列向量 转化为 列矩阵
# [[17.592  ]
#  [ 9.1302 ]
#  [13.662  ]
#  [11.854  ]
#  ......
#  [ 0.61705]]

theta = np.matrix(np.array([0, 0]))  # theta是一个(1,2)的矩阵，分别对应theta0和theta1
# [[0 0]]   转化为 1行2列 的矩阵
# 为什么都是0？因为我们从(0,0)点开始


# 看下维度
print(X.shape, theta.shape, y.shape)
# X：(97, 2)  theta：(1, 2)  y：(97, 1)


# 求代价函数J:
def computeCost(X, y, theta):
    # X 乘 theta 的 转置，减去y，再平方，得到方差
    inner = np.power(((X * theta.T) - y), 2)
    return np.sum(inner) / (2 * len(X))

print(computeCost(X, y, theta))
# 初次计算出的损失函数J的值为 32.072733877455676


# =====================================================================================================================
# 2、批量梯度下降
def grandientDescent(X, y, theta, alpha, iters):   # iters：迭代次数
    temp = np.matrix(np.zeros(theta.shape))  # 构造一行两列的0值矩阵
    parameters = int(theta.ravel().shape[1])  # ravel计算需要求解的参数个数，功能将多为数组降至一维 此处parameters=2
    cost = np.zeros(iters)  # 构建iters个0的数组

    for i in range(iters):
        error = (X * theta.T) - y
        for j in range(parameters):      #  j = 0 或 j = 1
            term = np.multiply(error, X[:, j])  # 计算两矩阵（h(x) - y）* x(j)
            temp[0, j] = theta[0, j] - ((alpha / len(X)) * np.sum(term))

        theta = temp
        cost[i] = computeCost(X, y, theta)  # 计算损失
        print('temp: ', temp, '  cost', cost[i])


    return theta, cost


# 初始化一些附加变量，学习速率α和要执行的迭代次数
alpha = 0.01  # 学习速率α
iters = 1000  # 1000次迭代


# 运行梯度下降算法，将我们的参数theta适合于训练集
g, cost = grandientDescent(X, y, theta, alpha, iters)
print(g)
# [[-3.24140214  1.1272942 ]] 即所求的theta0 和 theta1


# 再使用拟合的参数计算训练模型的代价函数（误差）
print("============================================")
print(computeCost(X, y, g))
#  4.515955503078912


# 现在来绘制线性模型以及数据，直观地看出他们的拟合，flg代表整个图像，ax代表实例
x = np.linspace(data.Population.min(), data.Population.max(), 100)  # 均匀抽100个样本
f = g[0, 0] + (g[0, 1] * x)  # g[0, 0]代表theta0， g[0, 1]代表theta1

fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(x, f, 'r', label="Prediction")
ax.scatter(data.Population, data.Profit, label='Traning Data')
ax.legend(loc=4) # 显示标签位置
ax.set_xlabel('Population')
ax.set_ylabel('Profit')
ax.set_title('Predicted Profit vs. Population Size')
plt.show()


#  绘制损失函数值的变化情况
fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(np.arange(iters), cost, 'r')  # 横坐标为0-1000，纵坐标为损失函数值
ax.set_xlabel('Iterations')
ax.set_ylabel('Cost')
ax.set_title('Error vs. Training Epoch')
plt.show()



# 3、多变量线性回归
data2 = pd.read_csv('ex1data2.txt', header=None, names=['Size', 'Bedrooms', 'Price'])
print(data2.head())
# 多变量的data2
#    Size  Bedrooms   Price
# 0  2104         3  399900
# 1  1600         3  329900
# 2  2400         3  369000
# 3  1416         2  232000
# 4  3000         4  539900

# 对于此任务，我们添加了另一个预处理的步骤 - 特征归一化
# 如果这个房子价格不归一化，它的数量级和你输入值归一化数量级差别太大，几十万的数量级和个位小数做回归，就不能保证收敛了，预测的y和实际上y几十万差的太多了
data2 = (data2 - data2.mean()) / data2.std()  # 每一维 减去均值，再除以方差
print(data2.head())
# 归一化后的数据如下
#        Size  Bedrooms     Price
# 0  0.130010 -0.223675  0.475747
# 1 -0.504190 -0.223675 -0.084074
# 2  0.502476 -0.223675  0.228626
# 3 -0.735723 -1.537767 -0.867025
# 4  1.257476  1.090417  1.595389

# 现在重复第一部分的预处理步骤，并对新数据集运行线性回归程序
data2.insert(0, 'Ones', 1)

cols = data2.shape[1]
X2 = data2.iloc[:, 0:cols-1]
y2 = data2.iloc[:, cols-1:cols]

X2 = np.matrix(X2.values)
y2 = np.matrix(y2.values)
theta2 = np.matrix(np.array([0, 0, 0]))  # 1*3 行矩阵

g2, cost2 = grandientDescent(X2, y2, theta2, alpha, iters)

print(computeCost(X2, y2, g2))
# 计算得到的损失函数的值为：
# 0.1307033696077189

# 查看一下这个训练过程
fig, ax = plt.subplots(figsize=(12, 8))
ax.plot(np.arange(iters), cost2, 'r')
ax.set_xlabel('Iterations')
ax.set_ylabel('Cost')
ax.set_title('Error vs. Training Epoch')
plt.show()


# 4、正规方程
def normalEqn(X, y):
    theta = np.linalg.inv(X.T@X)@X.T@y
    return theta

final_theta2 = normalEqn(X, y)
print(final_theta2)
